#!/bin/bash
set -e

while getopts  ":b:s:p:t:" opts
do
        case  $opts  in
		b)
			assembly=$OPTARG
			;;
        s)
			sample_name=$OPTARG
			;;
		p)
			out_prefix=$OPTARG
			;;
		t)
			total_reads=$OPTARG
			;;
		\?)
			echo `basename $0` '[-b assembly] [-s sample_name] [-p out_prefix] <in.bam> <([0.1|"0.1,0.3,0.5,..."])>'
			exit 65
			;;
		
        esac
done
shift $(($OPTIND - 1))


if [ $# -lt 1 ]; then
        echo `basename $0` '[-b assembly] [-s sample_name] [-p out_prefix] <in.bam> <([0.1|"0.1,0.3,0.5,..."])>'
        exit 65
fi

if [ -z "$assembly" ]; then
	export assembly=b37
else 
	export assembly
fi


. /mnt/ilustre/app/medical/tools/.var

reads_len=135
genome_len=3095677412
if test -z $total_reads
then
	if test -e $1.reads_count.txt; then
	total_reads=`cat $1.reads_count.txt`
	else
	total_reads=$(samtools view -c $1) #if input, save time
	echo $total_reads > $1.reads_count.txt
	fi
fi
echo; echo reads count: $total_reads

depth=$(echo ${total_reads} $reads_len $genome_len |awk '{printf("%.5f", $1*$2/$3)}')
echo;echo mean depth: $depth

fracs=$(echo $2 |cut -d',' -f1- --output-delimiter ' ')

for frac0 in $fracs; do
echo;echo frac: $frac0
frac=$(echo $frac0 $depth |awk '{printf("%.5f", $1/$2)}')
echo;echo frac: $frac
fract=$(echo $RANDOM $frac |awk '{printf("%.5f", $1+$2)}')
echo;echo fract: $fract

samtools view -@8 -F4 -b -s $fract $1 > $out_prefix.${frac0}x.bam
samtools index $out_prefix.${frac0}x.bam

summ_sc_gatk.sh -p $out_prefix.${frac0}x $out_prefix.${frac0}x.bam $data_path/intervals/1/cnv_sc/$assembly.rand.site.bed

done
echo done

exit 0